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Abstract 

We develop a new representation for the integrals associated with Feynman diagrams. This leads directly to 
a novel method for the numerical evaluation of these integrals, which avoids the use of Monte Carlo techniques. 
Our approach is based on based on the theory of generalized sine (sin(a;)/a;) functions, from which we derive an 
approximation to the propagator that is expressed as an infinite sum. When the propagators in the Feynman 
integrals are replaced with the approximate form all integrals over internal momenta and vertices are converted 
into Gaussians, which can be evaluated analytically. Performing the Gaussians yields a multi-dimensional infinite 
sum which approximates the corresponding Feynman integral. The difference between the exact result and this 
approximation is set by an adjustable parameter, and can be made arbitrarily small. We discuss the extraction 
of regularization independent quantities and demonstrate, both in theory and practice, that these sums can be 
evaluated quickly, even for third or fourth order diagrams. Lastly, we survey strategies for numerically evaluating 
the multi-dimensional sums. We illustrate the method with specific examples, including the the second order sunset 
diagram from quartic scalar field theory, and several higher-order diagrams. In this initial paper we focus upon 
scalar field theories in Euclidean spacetime, but expect that this approach can be generalized to fields with spin. 
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1 Introduction 

The expansion of relativistic quantum field theory devel- 
oped over 50 years ago by Feynman and Schwinger is still 
the basis for all perturbative calculations. However, an- 
alytical evaluations of graphs have become increasingly 
sophisticated, and digital computers make it possible to 
calculate graphs numerically. Since the accuracy of ex- 
periments in fundamental particle physics is continually 
improving, obtaining theoretical predictions whose preci- 
sion matches that of the data requires the evaluation of 
ever larger numbers of more complex graphs. Moreover, 
non-Abelian and effective field theories also increase the 
variety of diagrams that must be considered. Thus the 
evaluation of Feynman diagrams remains an active area 
of research. 

Many Feynman diagrams can be partially, or com- 
pletely, evaluated with purely analytical methods. When 
a diagram cannot be reduced to a closed form, the remain- 
ing integrals must be tackled numerically. Since these 
integrals may be multi-dimensional, Monte Carlo inte- 
gration is usually the only practical numerical algorithm 
0,^]. While conceptually simple, Monte Carlo methods 
converge slowly and accurate calculations of non-trivial 
graphs can be extremely time consuming. However, this 
combination of analytical and numerical techniques is ex- 
tremely powerful and, for instance, makes it possible to 



evaluate all the eighth order graphs in QED which con- 
tribute to the electron magnetic moment Q . In addition, 
computer algebra methods can be used to tackle the an- 
alytical categorization and simplification of diagrams, as 
summarized in the recent review Q . 

In this paper, our immediate goal is to present a new 
approach to the numerical evaluation of Feynman dia- 
grams. Its principal advantages are that it applies to ar- 
bitrary topologies, involves a minimal amount of analytic 
overhead, and has the potential to be much faster than 
Monte Carlo methods, especially when a result with more 
than one or two significant figures is required. The numer- 
ical method we develop is based on a novel representation 
of the Feynman integrals themselves, which allows us to 
approximate an arbitrary Feynman integral with a multi- 
dimensional infinite sum. 

The basis of our approach is an approximation to the 
spacetime propagator, derived using the theory of gener- 
alized Sine functions and expressed as an infinite sum. 
The approximation is governed by a small parameter, 
h, and converges rapidly to the exact result as h — > 0. 
For our purposes, the crucial feature of the approximate 
propagator is that its spatial (or momentum) dependence 
appears in terms like exp (— (x — y) 2 )- When the prop- 
agators inside a Feynman integral are replaced with this 
approximate form, integrals over vertex locations or inter- 
nal momenta required by the Feynman rules all reduce to 
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Gaussian integrals, which can all be performed explicitly 
for any diagram. Once the Gaussian integrals have been 
performed, the result is an iV-dimensional infinite sum, 
where N is the number of propagators in the diagram. 
We refer to this infinite sum as the Sine function repre- 
sentation of the Feynman integral. 

An immediate advantage of the Sine function represen- 
tation is that sums are much easier to compute numeri- 
cally than integrals, since sums are intrinsically discrete. 
More importantly, the general term in the sum decays 
exponentially, which greatly facilitates its numerical eval- 
uation. We show that the effort needed to evaluate the 
N dimensional Sine function representation of a Feyn- 
man integral increases linearly with the desired number 
of significant figures, although for a complicated topology 
the constant of proportionality can be large. This is in 
contrast to Monte Carlo methods where each successive 
significant figure is generally more costly to obtain than 
the last, even when the convergence is improved through 
the use of adaptive sampling. 

This paper introduces the Sine function representa- 
tion, and demonstrates its use in several explicit exam- 
ples. We also present an approach to regularizing di- 
vergent diagrams and discuss the extraction of renormal- 
ization independent quantities in the context of the Sine 
function representation. We focus our attention on scalar 
field theories in Euclidean spacetime, and will address the 
extension of Sine function methods to more complicated 
theories in future work. 

Finally, while the Sine function representation applies 
to perturbative quantum field theory, the approximate 
propagator is also a key ingredient of a new approach 
to "exact" numerical quantum field theory, the 'source 
Galerkin' method @-||. The source Galerkin method also 
eschews the use of Monte Carlo calculations. In addi- 
tion to improving its computational efficiency, by avoid- 
ing Monte Carlo techniques the source Galerkin method 
sidesteps many of the problems faced by other numeri- 
cal approaches to non-perturbative quantum field theory, 
particularly in models with fermions. While the physical 
bases of the Sine function representation and the source 
Galerkin method are very different, there is some intrigu- 
ing evidence that insights gained while constructing the 
Sine function representation may also allow us to improve 
the performance of the source Galerkin method. 

The structure of this paper is as follows: In Section 2 
we define the generalized sine functions, and review their 
relevant properties. The Sine function representation of 
the propagator is described in Section 3. In Section 4 
we show this form of the propagator leads to a new rep- 



resentation for Feynman integrals, and derive the "Sine 
function Feynman rules" for scalar field theory. We then 
apply these rules to the two loop "sunset" contribution to 
the \<f> two-point function, and compare our results to a 
conventional analytic calculation. In Section 5 we discuss 
techniques for efficiently evaluating the multi-dimensional 
sums obtained from the Sine function Feynman rules. We 
then evaluate the sums derived for representative third 
and fourth order diagrams, and compare these calcula- 
tions to direct integrations with VEGAS. Finally, in Sec- 
tion 6, we summarize our results and identify questions 
which we will pursue in future work. 

2 Generalized Sine Functions 

We begin by describing a generalized Sine function, which 
is defined by 



Sk(h,x) 



sin [n(x — kh)/h] 



(1) 



ir(x — kh)/h 

This is an obvious extension of the usual sine function, 
sinc(x) = sin(x)/x. Stenger ]l0| gives a thorough discus- 
sion of these functions, while a more introductory account 
is provided by Higgins jll] . We follow Stenger in using the 
capitalized "Sine" to distinguish the generalized version, 
although our notation for Sk{h, x) differs from his.[] 
The Sine function has the integral representation, 

S k (h,x) = ^- k e^-^dt, (2) 



2?r 



-ir/h 



which is trivial to prove. We recognize equation (Q) as 
the Fourier transform of a finite wavetrain, which is non- 
zero only in the interval [—Tt/h,Tt/h]. Using the integral 
representation, we can establish the orthonormality of the 
Sine functions, 



dxSk(h,x) = h, (3) 
dx Sk(h,x)Si(h,x) = hdki- (4) 



Any function f(z), which is analytic on a rectangular 
strip of the complex plane, that is centered on the real 
axis with width 2d, is approximated by 



/(*)« J2 f(kh)S k (h,z). 



(5) 



1 Our Sk(h,x) is equivalent to Stenger's S(k, h)(x). 
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The approximation improves as h decreases, and this re- 
lationship is quantified below. Equation ([5]) is referred to 
as the "Sine expansion" or "Sine approximation" of f(z). 
The Sine expansion is exact for Paley-Weiner, or band- 
limited, functions, which corresponds to f(z) be repre- 
sentable in the following way: 



f(z)= T g(uy zu du, geL 2 (-7r/h,Tr/h). (6) 
J-w/h 



This states that / is a function whose Fourier transform, 
g, has compact support on the interval [—ir/h, ir/h] or, al- 
ternatively, that / has no frequency outside the "band" , 
[—Tr/h,ir/h]. Sine functions are used in a wide variety of 
areas, and Stenger gives a number of examples from ap- 
plied mathematics and classical physics. They are also 
at the basis of Shannon's Sampling Theorem and are fre- 
quently encountered in communications theory. 

For our purposes the most important property of the 
Sine function is that if / is approximated by its Sine ex- 
pansion, equation (|^), we can quickly derive a related ap- 
proximation for the definite integral, 



f(z) dz 



00 r oo 



I dzf(kh)S k (h,z) 

00 

h J2 f(kh) (7) 



k— — 00 

where the last equality follows from the normalization, 
equation (||) . In the following section we use this relation- 
ship to approximate the propagator as an infinite sum. 

From Theorem 3.1.2 of Stenger, the error in the Sine 
approximation to a function f(z) which is analytic in an 
infinite strip of the complex plane, centered upon the real 
axis with width 2d, is 



5f(h,z) 



f(z)- J2 f(kh)S k (h,z) (8) 
f(t - id) 



k=—oo 



sin 2^ roo 
h 



2-ki 



(t-z- id) sin [f(t - id)] 
f(t + id) 



(t-z + id) sin[f(t + ZG0] 



dt.(9) 



The discrepancy between J"f 
tion is 



/ and its Sine approxima- 



A/(fc) 



dzf(z) 



OO 

' E 

k— — oc 



f(kh). 



(10) 



Clearly, Af(h) is the result of integrating Sf(h, z) be- 
tween ±00, and Theorem (3.2.1) of Stenger allows us to 
deduce that 

' A ^'^2^h^ 
where C is a real number which is independent of h. When 
h < d we see that |A/(/i)| drops exponentially as h de- 
creases linearly. 

3 The Scalar Field Propagator 

3.1 The Spacetime Propagator 

Feynman diagrams are conventionally computed in mo- 
mentum space, since the propagators are algebraically 
simpler, and the momentum conservation constraint at 
each vertex reduces the number of integrals that would 
otherwise need to be performed. However, we will find it 
more convenient to perform most of our work in coordi- 
nate space. 

The spacetime propagator for a scalar field theory is 
the Fourier transform of the momentum space propagator, 



namely^ 



G(p) 



G{x -y) = 



1 



d%p 



a ip(x-y) 



(2tt) 4 p 2 + TO 2 



(12) 



(13) 



where x, y and p are Lorentz four-vectorsjj and we are 
working in Euclidean space. Redefining x — y as x and 
adding a cut-off in anticipation of divergent diagrams, we 
can write 



G A (x) 



d%P 

(2tt) 4 p 2 



V/A 2 



Exponentiating the denominator gives 

1 

G A (x)= I ^rre* px -P^ A ~ I dse~ s ^ +m2 \ 



(2k) 



(14) 



(15) 



A rationale for this regulated form was given by Hahn [g| ; 
in particular, formulating the cut-off in this way pre- 
serves the Gaussian form of the integrals (as would a 



2 This statement also implicitly defines the normalization we have 
adopted for the Fourier transform and its inverse. We also assume 
m 2 > 0. 

3 We will suppress the Lorentz indices on four- vectors unless they 
are directly relevant to the calculation at hand. 
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dimensionally- regulated formulation). Reversing the or- 
der of integration, completing the square, performing the 
Gaussian integrals, and finally shifting from s to s/m 2 
yields 



G A (x) = 



ds ■ 



exp 



4(s+fr) 



(16) 

Making a further change of variable to, s = e 2 , we obtain 

7T1 2 

G A (x) = — — — / dz 



(4tt) 2 



x exp 



4(e 2 + fr) 



(17) 



The crucial step is to approximate G A (x) as an integral 
over a Sine expansion Q of the form given by equation (]?]) , 



G Ah (x) 



(At:) 2 



/+oo jx> 
-OO j. 



„kh 



k= — oo 



x exp 



„kh 



(e kh + fr) 2 

^,2^,2 



4{e kh + %) 



S k (h,z) 



i 2 h 



(4tt) 2 



E 



,-,kh 



'pkh i \2 



x exp 



-e kh - 



4(e 



A 5 ' 



(18) 



When A — > oo the coordinate space propagator, G(x), can 
be expressed as a modified Bessel function of the second 
kind, 

m 2 -fTi(w|:r|) 



G(x) 



(19) 



47r 2 m|x| 

In this limit, the Sine approximation to G A also simplifies, 
m 2 h 



G h (x) 



(4tt) 2 



E ex P 



fe=-c 



-kh — e 



(20) 



Finally, it will be convenient to introduce a more concise 
representation for the various factors that appear in the 
general term of sum, 



G Ah (a - b) 
where 

c(fc) 
p(k) 



mrh 
jAnj 2 



^2 p(fc)exp 



m 2 (a - 6) 2 
4c(fe) 



= e feh + 



TTi 

A 2 " 



e fe "- exp (— e ) exp (fc/i — e fc/l ) 



(e** + fr) 2 



c(fc) 2 



(21) 

(22) 
(23) 






01 




-10 


-L u 






-18 


10 






-2 6 


10 






-34 


10 






-42 


10 





0.2 



. 4 



. 6 



. 8 



1. 



Figure 1: The relative error of the approximation G^ is 
shown, by plotting \AGh/G\ as a function of h for m = 
x = 1 and A — > oo. 

We now have four different versions of the propagator, 
which are related as follows. Firstly, the exact propagator, 
G(x), is given by equation (|l9|), while equation ( |l6| ) de- 
fines the cut-off propagator, G A (x). These expressions are 
approximated by Gu(x) and G A h{x) respectively, which 
are defined by equations ( p(i| ) and (|l^). As we shall see, 
the approximations rapidly approach the exact values as 
h becomes small. 

Both the cut-off propagator, G A (a — b), and its ap- 
proximation, G\h(o- — b), depend only on the combination 
(a — b) 2 , and thus neither the introduction of A nor the 
Sine function approximation has violated Poincare invari- 
ance. Moreover, by Fourier transforming G A h{x) we could 
obtain a Sine function approximation for the momentum- 
space propagator, equation (12), in the presence of the 
cut-off. In Section 4, we develop the Sine function Feyn- 
man rules in coordinate space, but the momentum space 
rules could be obtained using an almost identical argu- 
ment. 

3.2 Accuracy of the Approximation 

Our next task is to assess the accuracy of the approxima- 
tion involved in writing Gh{x) and G A h{x). We define 



AG Ah (x) = G A (x) - G Ah (x). 



(24) 



The Sine function approximation to the propagator is de- 
rived from equation (0) , so the h dependence of the error 
is governed by equation (|ll]). The form of f(z) which 
appears in the left hand side of equation ([ll]) is our inte- 
grand, 



/(*) 



■ exp 



4(e^ + ^) 



(25) 



4 



N 

-9 
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Figure 2: The relative error of the approximation Gh is 
shown, by plotting \AGh/G\ as a function of x for m = 1, 
/t, = 0.4 and A — > oo. 

where x, m and A are simply parameters. Since f{z) is 
infinite when e z +m 2 /A 2 =0, or when z = 2 log(r7i/A)±i7r, 
it follows that d, the width of the strip of the complex 
plane which appears in equation (^), can take any value 
in the open interval (0,7r). Using the maximal possible 
value of d, we derive the following bound on the accuracy 
of the approximation, 

\AG Ah \~e- 2 */ h , h^O. (26) 

Thus as h approaches zero, Af(h, z) decreases expo- 
nentially. Looking at the general form of the Sine expan- 
sion, we see that the magnitude of the general term in 
the sum decreases (at least) exponentially with increas- 
ing k. Since any numerical evaluation of these sums must 
be cut off at finite |fc|, it follows that while G^h will ap- 
proach G\ exponentially quickly as h goes to zero, the 
number of terms that need be considered when comput- 
ing the sum grows (at worst) linearly with 1/h. This im- 
pressive convergence facilitates the numerical evaluation 
of the approximate propagators and, consequently, the 
Sine function representations of the Feynman integrals. 

Having discussed the abstract convergence properties 
of the Sine approximation to the scalar field propagator, 
we now present some specific numerical results. Fig (Q) 
shows AGfc(l)/G(l) as a function of h. These results 
were obtained using arbitrary precision arithmetic. When 
h = 0.1, the two expressions differ by less than 1 part in 
10~ 40 , and if h is even slightly less than unity the accu- 
racy of the approximation is phenomenally impressive. As 
a specific example, choosing h — 0.5 gives an accuracy of 
1 part in 10 8 , while h = 0.25 yields an agreement better 



Figure 3: The number of significant terms in the series for 
Gh(l), with m = 1 is displayed. When the double infinite 
series for Gh(x) is truncated at k — —N and k — N, the 
result is accurate to within approximately 1 part in 10 14 . 
between 

than 1 part in 10 16 . On most computers, double preci- 
sion real numbers contain 16 decimal digits, so if h < 0.25, 
Gh{x) and G(x) will be numerically indistinguishable. In 
fig. (||) we show that for fixed h, the accuracy of the ap- 
proximation is essentially independent of x. 

In figs (Q) and (|J) we have taken A — > oo for conve- 
nience, since this allows us to use the exact expression 
for the propagator, equation (|l9|). However, using a fi- 
nite value of A makes no significant difference to these 
conclusions. 

Our next concern is to determine how many terms in 
the series representation of Gh contribute significantly to 
the final result when it is evaluated numerically. Fig (|^) 
shows the number of terms in the series for Gh which 
contribute more than 1 part in 10 14 of the final result. If 
h = 0.25, we need to include terms for which |fc| < 21. 
Note that this measures the convergence to Gh rather 
than G, but for h < 0.25 the difference is less than the 
truncation error in the numerical calculation of Gh{x). 
Also, if we do not need a highly accurate calculation, we 
can increase the value of h and thus reduce the number 
of terms that contribute to the sum. 

4 Evaluating Feynman Diagrams 

4.1 The Sine Function Representation 

Having developed the Sine function approximation to the 
propagator, we now turn to the evaluation of Feynman 



5 



diagrams. In coordinate space the Feynman rules specify 
that we construct the product of the propagators corre- 
sponding to the N internal lines of the diagram, and then 
integrate over each of the M internal vertices. If the ex- 
ternal lines are attached at the points X\ , X2, ■ •• , %a an d 
the internal vertices are labeled by y±, y%, • • • , yu we are 
therefore integrating over all spacetime for each of the yj . 
Up to overall prefactors consisting of powers of the cou- 
pling constants and the topological weight of the diagram, 
these integrals take the form 



Ia 



d 4 yx ■ ■ ■ d 4 y M G A ( Xl - yi) ■ • • 

G\{x A - yj) ■ ■ ■ G A (y k - vm)- (27) 



In general, the I A cannot be reduced to a closed, analytic 
form. We now use the Sine expanded propagator, G A h, to 
derive the Sine function representation, I A h , of a Feynman 
integral, I A . 

Firstly, G A (x) — G A h{x) + AGa/,(i). For given A and 
h, there exists a number c, such that AG A h(x) < cG A (x), 
and c <C 1 if h is sufficiently small. We therefore obtain 

Ia = J d 4 2/i • • • d 4 y M [ {G Ah (xi ~ 2/1 ) + 

AG A i(ii - y\)) ■ ■ ■ 
[GAhiVk - Vm) + AG Ah (y k ~ yu))\ (28) 

d 4 yi ■ ■ ■ d 4 y M G A h(xi - J/i) • • • G Ah (y k - Vm) 

+ M J d 4 yi ■ ■ ■ d 4 y M [AG Ah ( Xl - Vl ) ■ ■ ■ 

G Ah {y k -y M )], (29) 

where we have assumed that h is small enough to make 
AGaa(i) «C G A h(x), so we can drop terms that are sec- 
ond order in the AG A h(x). Consequently, we can use the 
definition of c to write 



Ia 



d 4 yi ■ ■ ■ d 4 y M G Ah (xi -y%)--- G Ah (y k - Vm) 

(30) 

(31) 



+cM J d 4 yx ■ ■ ■ d 4 y M G Ah (xi - yi) ■ 
G Ah {yk - vm) 

lAh 



where 

lAh = J d 4 yi ■ ■ ■ d 4 y M G Ah {xi - yi) ■ ■ ■ G Ah (y k - vm)- 

(32) 

We see that the Sine function representation, I A h, dif- 
fers from I A by 0(c). Since c — > as h — * 0, I A h is an 
arbitrarily accurate approximation to I A . 



The spatial dependence of G A h(yi — y-i) occurs only 
in terms like exp (— (yi — 2/2) 2 )- Thus the integrals over 
the internal vertices in equation ( |32"|) are all reduced to 
Gaussians, which can be performed analytically, no mat- 
ter how complicated the diagram. Thus the Sine function 
representation is an N dimensional infinite sum, where N 
is the number of internal lines in the diagram. In general, 
this sum has the form 



I Ah = h N y^J(kx, ■ ■ ■ ,k N ,Xi,- ■ ■ ,x M , h), 



(33) 



where the stands for the N individual summations 
between ±00, corresponding to each of the k\. Referring 
back to equation (^l|), we see that the fcj only appear in 
/ inside the c(fci) and p(ki) of the original propagators. 
From the form of G A h and the properties of Gaussian 
integrals, it follows that that the spatial dependence of 
/ is restricted to terms with the form exp (—(x a — xi,) 2 ). 
Consequently, like G A h{x) the I A h(x) are also manifestly 
covariant. 

For reference, a Gaussian integral in four (Euclidean) 
spacetime dimensions has the general form 



d 4 x exp (—a + 2b^x^ — dx^) 



d 2 



■ exp 



(34) 

and we have actually already used this result to derive 
G A (x). Furthermore, the Fourier transform of a Gaussian, 



2 / 2 

n ( p 



d 4 xe~ lpx e- dx = -^exp . 

cr \ Ad 



(35) 



is simply a special case of equation (|34|). Using this result, 
it is easy to transform I A h (x) , and obtain the correspond- 
ing momentum space representation. 

Before we discuss the convergence properties and eval- 
uation of these sums, we note that the spatial integrations 
in I A can also be performed if we work with integral rep- 
resentation of the cut-off propagator, equation (|l6|). Con- 



sider equation (27) 



I A = 



m \ 
4W 



2N 



Si ■ ■ ■ S N 



/^..^ex P (~^ 



N -a- 

n e 

m 2 (x 1 - yi) 2 \ 
+ £) ) 



exp 



m 2 (y k - y M f 
' 4( SW + ^) 



(36) 



Combining the Gaussians, we obtain an TV-dimensional 
integral over the parameters, s,-. It is this integral which 



G 



is directly approximated by the Sine function representa- 
tion, Iah- 

Writing an iV-dimensional integral as an N dimen- 
sional sum is not necessarily an improvement, since the 
sums still have to be numerically evaluated. Moreover, 
while I\h is an iV-dimensional sum, other approaches ex- 
press to the Feynman integrals with less than N dimen- 
sions. For instance, if we had integrated over the exact 
form of the scalar propagator in 4 — e dimensions the re- 
sult would still be finite, but we only need to perform 
M integrals over the internal vertices. Obviously, M is 
always less than N, and the Euclidean integrals can typ- 
ically each be expressed as a single integral over a radial 
variable, even though they are not Gaussians. Likewise, 
in momentum space the dimension of the integral is de- 
termined by the number of independent momenta, which 
must be also be smaller than the number of propagators. 
However, it will turn out that the Sine approximations 
to Feynman integrals inherit the rapid convergence prop- 
erties possessed by Gm,, and even complicated diagrams 
are computationally tractable. 

At this point it is useful to consider a naive estimate of 
the computational cost of evaluating I a using the infinite 
sum derived from I ah- If the diagram contains N propa- 
gators, and we truncate the individual sums at |fcj| ~ n, 
we would expect to evaluate ~ n N terms. A complicated 
diagram can have N ~ 10, and in this case modest values 
of n will still make n N exorbitantly high, even when we 
are evaluating millions of terms a second. However, in the 
next sections we show that this simple calculation is actu- 
ally far too pessimistic. Firstly, the number of significant 
terms grows more slowly with N than the volume of an 
A-dimensional hypercube, and a sensible numerical algo- 
rithm will focus on these terms. More radically, we also 
consider approaches which reduce the overall "power" of 
the problem to a number less than TV. 

Even at this point, though, we can probe the relation- 
ship between the accuracy of the approximation to I a and 
the number of terms that must be evaluated in I Ah ■ Re- 
call that AG Ah decreases exponentially with ft, -1 , while 
the number of terms significant terms in the sum increases 
linearly. Since the number c that appears in the deriva- 
tion of equation is effectively the maximum value of 
AGAh(x) for < x 2 < oo, I ah inherits the convergence 
properties of Gaii(x). Consequently, to double the num- 
ber of significant digits of I a given by I Ah > we must halve 
the value of h. Roughly speaking, this increases the num- 
ber of terms which contribute to I Ah by a factor of ~ 2 . 
This 2 N may seem a little daunting, but the important 
point is that the number of significant terms still grows 



linearly with the desired accuracy. 

It is this relationship between h and the number of 
significant terms in lAh that holds out the promise of a 
dramatic improvement over Monte Carlo methods. For 
comparison, while adaptive Monte Carlo routines such as 
VEGAS focus on the most important parts of the over- 
all volume, obtaining high accuracy is inherently difficult 
because of the statistical nature of the method. At worst, 
the accuracy of the Monte Carlo algorithm scales as 1/ y/ri 
where n is the number of points to be evaluated. Since 
each decimal place in the result increases the accuracy by 
a factor of 10, the CPU time needed for each successive 
significant figure can be as much as 100 times greater than 
that required for its predecessor. Thus the time required 
by a Monte Carlo integration can rise exponentially with 
the desired accuracy, in contrast with the linear increase 
that applies to the Sine function method. 

4.2 Renormalization and Sine 
Function Methods 

The next issue we discuss is regularization and renor- 
malization when using the Sine function representation. 
Firstly, we wish to remove a possible source of confu- 
sion by clarifying the two different limits implied by Iahj 
namely A — > oo and h — > 0. In the limit h — ► 0, the error 
implicit in the Sine function representation drops to zero. 
However, when we use the Sine function representation 
as the basis of a numerical calculation, we do not (and 
cannot) evaluate this limit exactly. Rather, we choose a 
finite value of h that is small enough to ensure that the 
error induced by the Sine expansion is less than the accu- 
racy we require from our computation. Beyond this, the 
value of h has no physical significance. 

In contrast, when A — > oo most non-trivial integrals 
diverge, and the corresponding infinities must be removed. 
Also, recall that while Feynman integrals can be regular- 
ized (rendered finite) in many different ways, renormal- 
ized quantities cannot contain any residual dependence 
on the regularization scheme. 

In this initial paper we focus on diagrams which con- 
tribute to the propagator. The possible divergences of a 
diagram with no sub-divergences scale as p° and p 2 , so for 
a given £(p) which contributes to the propagator we can 
form a renormalized quantity, 



E(p) = lim 

A^oo 



E(p)-£(0)-p 



2 ^(P) 



d{p 2 ) 



p=0 



(37) 



where the two terms that have been dropped are sim- 
ply the first two contributions to the Taylor expansion 
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of the regularized diagram, S(p 2 ). We can see no a pri- 
ori reason not to use dimensional regularization, but in 
practice we have found the covariant, A-dependent cut- 
off to be convenient. However, any cut-off scheme which 
did not ensure that the integrals over the internal vertices 
are always Gaussians would undermine the Sine function 
representation, since there would be no guarantee that 
the spatial integrals in an arbitrary diagram could all be 
performed. 

A given S(p) which contributes to the propagator has 
a Sine function representation with the generic form 



S(p)~^Aexp(-/5) 



(38) 



where A and B are functions of the fcj, A, h, and to, but in- 
dependent of p. In principle, we could evaluate the corre- 
sponding by computing the terms on the right hand 
side of equation (|37|) at a large (but still finite) value of 
A, and then performing the subtractions. Unfortunately, 
this naive approach is foolish since the individual terms 
in equation (37) are dominated by the divergences and 
combining them would lead to a drastic loss of numerical 
precision. 

We resolve this problem by writing the right hand side 
of equation (|37|) as a single multi-dimensional sum, 



~ V I" lim A{exp(-p 2 B) - 1 + p 2 B) 

L — ' LA— >oo 



(39) 



Mathematically this is not objectionable, since by combin- 
ing the individual terms while keeping A finite we never 
manipulate formally divergent sums. We take the limit 
A — > oo in equation ( |39"| ) analytically, before we compute 
the multi-dimensional sum. Not only does this avoid the 
loss of precision inherent in the naive evaluation of equa- 
tion ([37|) , it yields an analytical, renormalized, cut-off in- 
dependent expression for S(p). 

The numerical task we face is to evaluate the sum that 
results from taking the limit in equation (|39|). In general, 
this is no more difficult than evaluating the cut-off depen- 
dent expression, equation (38). The one caveat is that if 
p 2 B< 1 the individual terms in the summand of equa- 
tion ([39]) have very similar magnitudes, and combining 
them directly would lead to a loss of numerical precision. 
However, when p 2 B <C 1 the combination we need can be 
computed accurately by replacing the exponential with 
the first few terms of its Taylor series, which converges 
very quickly. 

We give a specific example which illustrates this ap- 
proach to renormalization when we compute the sunset di- 
agram from Xcj) 4 . Moreover, we believe that this argument 



can be generalized to diagrams with sub-divergences, but 
we will reserve a detailed discussion of this problem for a 
future publication. 

4.3 Sine Function Feynman Rules 

We are now in a position to write down the analog of the 
conventional Feynman rules that specify how to obtain 
the Sine function representation for an arbitrary diagram 
in scalar field theory: 

1. Write down the integral to be evaluated, using the 
usual coordinate space Feynman rules and propaga- 
tors. 

2. Replace the propagators with G\h(x), the Sine ex- 
pansion of the cut-off propagator. 

3. The spatial integrals are now reduced to Gaussians, 
which are performed analytically. 

4. Once the Gaussian integrals have been performed, 
Fourier transform the result into momentum space 
(if desired). 

5. Extract the regularization independent quantities 
(if desired). 

6. Evaluate the resulting multi-dimensional sum nu- 
merically. 

Before proceeding, we make two observations. Firstly, the 
Gaussian integrals can easily be performed using com- 
puter algebra packages, and the result expressed in the 
syntax of compiled languages such as Fortran or C. This 
suggests that the derivation and evaluation of the Sine 
function representation for a given group of diagrams will 
be easy to automate. Secondly, while the general term 
in the sum will grow more complicated as the number of 
internal vertices is increased, the integrals that must be 
performed are always Gaussians. Consequently, extract- 
ing the Sine function representation is not intrinsically 
more difficult for higher order diagrams than it is for sim- 
ple diagrams. 

4.3.1 Case Study I: The Sunset Diagram 

To give a specific example of our method, we evaluate the 
sunset diagram which contributes to the propagator of 
A0 4 scalar field theory, and which is depicted in fig. (Q). 
This diagram is a useful test subject, since its overlap- 
ping divergences ensure that it is has non-trivial structure, 
but it is still simple to enough to be treated analytically. 
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Figure 4: The sunset diagram contribution to the A</> 4 
propagator, with momentum, p. 



Consequently, we can use it to to test the Sine function 
method to an arbitrary level of precision. 

In coordinate space the diagram has a particularly 
simple form, since there are no internal vertices to in- 
tegrate across. Consequently, we have 



£ Aft (as - V) = G Ah (x - y) 3 



(40) 



which is just the product of three propagators. We use the 
same convention here as we do with G A h and Ixhi where 
the cut-off dependent Sine approximation is denoted by 
the subscripts A and h. We have suppressed the prefactor 
consisting of A 2 and the topological weight, since it plays 
no part in the analysis. Working through the prescription 
given above, 




Figure 5: The contribution of the sunset diagram to a is 
plotted. The data shown here was obtained from the Sine 
approximation with h = 0.4, but is indistinguishable from 
the exact result. The initial factor of m 2 /(47r) 4 in E(p) 
has been omitted from these results. 



T, Ah (x~y) = 



exp 



m 2 {x — y) 



v c(fci) c(k 2 ) c(k 3 
After carrying out the Fourier transform, we obtain 

p(k 1 )p(k 2 )p(k 3 ) 



(41) 



^Ah(p) 



i 

^77 



exp 



1 



c(fc 3 



i 



(42) 



We have now turned the integral represented by the sunset 
diagram into an infinite sum over all possible {^1,^2,^3}. 
Starting from equations ( |l6| ) and (|36|), we can derive the 
triple integral to which Ea/i is the approximation, 



Sa(p) 



(4tt) 4 



dsids2<ls3 Y\ 
1 



1=1 



(sT 1 - 

exp 



2 — 1 



*2 + So 



(43) 



where Si stands for Si + m 2 / A 2 . The parallel between Ea?j 
and Sa is clear. 

A little analysis establishes that both Ea and Ea/j 
are finite, but that they diverge as A — * 00, which is 
precisely what we expect, since this diagram contains an 
(overlapping) divergence. Analyzing the general term in 
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Figure 6: The relative error in the Sine function calcula- 
tion of the finite part of the sunset diagram, S, is shown, 
compared to obtained from the analytical result, equa- 
tion (H), with h = 0.4. 
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Figure 7: The relative error in the Sine function repre- 
sentation of the sunset diagrams is plotted as a function 
of h, with a fixed p of 1.4. Even when h ~ 1, the Sine 
function method returns several significant figures of the 
exact result, while the nested sum can be evaluated in a 
fraction of a second. 



the sum, we see that the divergence is "worst" when we 
sum along the line k\ — ki — k$, k\ < 0. Applying the 
renormalization prescription given by equation ( |39| ) we 
form 



771 

(4 



'•h 3 



— Y : 

-) 4 V fi 



p(kx)p(k 2 )p{k 3 ) 



exp 



2 



P 



i 



i 



c(fci) ^ c(fc 2 ) + c(fc 3 ) / 

The result of evaluating equation is plotted in fig. (||). 

We now compare this result to a standard analytic 
computation. Rather than tackle the evaluation of £a 
directly, we use the results of Ramond |l2j , who works in 
4 — 2e dimensions to compute that 



E(p) = 



-1 
1 - 2e 



(» 2 ) 2e (3m 2 K( P )+ P »K ll (p)) (45) 



where \x is the renormalization scale and the K(p) and 
p^K^lp) can be expanded as functions of e to give 



K(p) = 



r(2e) 1 



(4tt) 



4-2e 



l + e(l-21ogm 2 ) + 



7T" 

g" 

2 j'dx f 1 dy(l-y)\og{y) d]inl>i] 
lo Jo 

-2 / da; / dj/log(/)log(j/) 



e 2 ( 3 - + 21og 2 (m 2 ) - 41og(m 4 



r/i/ 



o 







(46) 



where 



/ = p 2 y( 1 - v) + m2 1 - V 



x(l — x) 



(47) 



In the two unevaluated integrals, a single integration can 
be performed without too much difficulty, and the result 
expressed in terms of dilogarithms. In addition, 



P^K^p) = 



r(2c) i 



(4tt 



i4-2e 



1-eQ+log (to 2 ) 



/?+ 1 



where 



-2/ d^l-^lo^^ 



TO 2 y 

P 2 ( 1 -y)?/ + m2 ( 1 -2/)' 



(48) 
(49) 
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Normally, infinities are eliminated from dimensionally 
regularized diagrams by dropping terms proportional to 
1/e 2 and 1/e. However, if we drop the first two terms in 
the Taylor expansion with respect to p 2 , and then take 
e — > we remove the divergent terms, and also ensure 
that the result is (in theory) equivalent to equation (44), 
in the limit h — > 0. The agreement between the exact 
result and the Sine function approximation is shown in 
figs ([j]) and (Q) . As predicted, we can obtain an arbitrarily 
good agreement with the exact result by lowering h, but 
we also find that discrepancy between the approximation 
and the exact result remains small even when h is of order 
unity. Larger values of h reduce the number of terms that 
contribute significantly to the sum, and if we only want a 
few significant figures we can set h *~ 1, which makes the 
diagrams very easy to evaluate numerically. 

5 Numerical Implementation 

5.1 A Simple Approach 

So far, we have focused on the theoretical development 
of the Sine function representation. In this section we 
turn our attention to the practical issues that arise when 
numerically evaluating the sums, and discuss the CPU 
time needed for typical calculations. 

Naively, computing a sum is a much simpler process 
than evaluating an integral since a sum is inherently dis- 
crete, whereas an integral is a continuous object which 
must be discretized before being numerically evaluated. 
The one proviso is that we must truncate the sum care- 
fully, to ensure that terms which do not contribute signifi- 
cantly to the final result are not needlessly calculated, and 
that terms which do not contribute are not accidentally 
ignored. 

In general, we are faced with the task of computing 
an A-dimensional infinite sum. For the case of diagrams 
which contribute to the propagator, the general term has 
the form of equation (|38|), and the overall sum can be 
written 



E 

fei 



few 



, (50) 



where / implicitly depends on m, p, A and h. We have 
focused on the renormalized sunset diagram, but the fol- 
lowing discussion is also applicable to the case with finite 
A. When computing the sums which represent Feynman 




Figure 8: Terms which contribute to contribute to S(p), 
with h = 0.4 and m = p = 1, by more than 1 part in 
10 10 are plotted. To evaluate the sum efficiently, we must 
focus on {fci, ks} inside this volume. 



integrals, it may not be practical to determine a priori 
which {fci, ■ • • , fejv} are significant. Consequently, we have 
developed a simple algorithm that evaluates I\h by dy- 
namically tracking the size of the general term as a pro- 
portion of the "running total" for each of the N nested 
sums in the overall sum. The precise criteria used to 
select the significant terms must be determined heuristi- 
cally. The truncation error must be watched carefully, as 
all the terms in lAh are positive, and truncating the sum 
always underestimates Ihh{p)- The volume of {fci, fe, £3} 
space which contributes to is shown in fig. (|J). 

The ki only appear in /(fci, • • • , k^) via the c(fcj) and 
p(ki), defined in equations ( |2^ ) and (23). Because m 2 
(and A 2 ) does not change during the course of a given 
calculation, the c(fc) and p(k) depend only on the integer 
values of the individual fcj. Thus we compute c(k) and 
p(k) for a range of ki and storing these results in arrays 
before starting on the evaluation of the sum itself. By 
using these stored values when the calculating the gen- 
eral term of the sum we avoid the need to recompute the 
exponentials in the p(k) and c(k) at every step. 

We start each successive sum with an initial value of fej 
which maximized the general term for the previous sum. 
This is efficient since only differs by ±1 between the 
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two cases, and the dependence of the general term on any 
of the ki is sufficiently weak to ensure that we are starting 
close to the maximal value. The one drawback to this ap- 
proach is that when computing the final terms in the sum, 
we are adding small numbers to a large one (the "running 
total" for the overall sum), a circumstance that can lead 
to a loss of accuracy in fixed precision arithmetic. We 
have ensured we are not losing accuracy in this way by 
spot-checking our results with quadruple-precision arith- 
metic, but in a more sophisticated code this issue would 
be addressed more elegantly. 

In the case of the sunset diagram, evaluating equa- 
tion (Q) with h = 0.4 takes around 3 seconds of work- 
station CPU time, and the final result differs from the 
exact value determined from equation ( f45|) by 1 part in 
10 10 . For h = 0.6, the accuracy drops to a few parts in 
10 s , and the execution time drops by around a factor of 
3.0 For the specific case of the sunset diagram, consider- 
able progress can be made with the analytic evaluation of 
the Feynman integral and this knowledge can be used to 
produce numerical values for However, when com- 

pared to other methods for numerically evaluating inte- 
grals, the convergence properties and speed of the Sine 
function method is impressive, involving as it does the 
direct numerical evaluation of a triple integral. 

5.2 Improvements to Convergence 

While simply adding up all the terms larger than a given 
minimum size is the most obvious way to evaluate an 
infinite sum, it is not necessarily the most efficient. In 
this initial paper we will give only a brief survey of tech- 
niques that could improve the convergence of the sums. 
In general, we foresee three broad classes of technique for 
improving convergence; analytical manipulations of the 
infinite sum's general term (e.g. pjj), numerical extrapo- 
lations of successive partial sums (e.g. [|[^4)), and meth- 
ods that effectively reduce the order of the N dimensional 
sums generated by the Sine function Feynman rules. 

We will not attempt analytical rearrangements or eval- 
uations of the sums we have derived. However, an inter- 
esting line of enquiry would be to consider whether known 
analytical treatments of Feynman integrals have analo- 
gous results which can be used to evaluate at least some of 
the infinite series in the I\h- We do consider one example 
of the other two techniques - namely using extrapolation 



4 More detailed information about the numerical calculations is 
given in an appendix. 
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Figure 9: The general term of the Sine integral form of 
the sunset diagram is shown, for h = 0.4, m = p = 1 
and ki = k^ — k% = k. The general term decays like a 
power-law as k becomes large and negative. 

to speed the evaluation of the innermost sums, and gen- 
crating a "dictionary" of pre-computed values from which 
the innermost sum can be accurately interpolated. 



5.2.1 Aitken Extrapolation 

For a general /atm if any one of the ki is becomes large 
and positive, the general term is dominated by factor like 
exp(— e kih ). This expression decreases so precipitously 
that the transition between terms which are significant 
and those which are negligible is very sharp, and there is 
little to be gained from extrapolative techniques. How- 
ever, for ki <C 0, the general term takes on a power- 
law behavior, whose regularity can be exploited by ex- 
trapolation. The general term of the sum is shown (for 
h = k 2 = h) in fig. @. 

As a specific example of numerical extrapolation meth- 
ods, consider the Aitken S 2 approach M, which is appli- 
cable to infinite series whose general term behaves like 
~ x k for large k, which is true of our problem when ki 
is decreasing. For a general infinite series with three suc- 
cessive partial sums S n —i,S n ,S n +t, Aitken's S 2 method 
gives the improved estimate, 



S — S n 



(S„ 



S n ) 2 



S n +1 2Sn + S n — 1 



(51) 



If the general term of the series is f(n), we can rewrite 
this as 

S' = S n+1 + { {n+ f( 1] " v (52) 
f(n) - f(n+l) 

Written in this form, we can recognize the general result 
for a sum of the type X^n 2 -™' with f(n) — x n . For 
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the specific case of the sunset diagram, using the Aitken 
S 2 extrapolation as k 3 becomes large and negative can 
improve the performance of the algorithm by at least 25%. 

5.2.2 Reduction of Order 

As we noted in the previous section, the Sine function 
Fcynman rules yield an A-dimensional infinite sum, where 
N is the number of internal lines in the diagram. Conse- 
quently, if we can reduce the dimension of the infinite sum, 
we can expect a significant improvement in efficiency. By 
considering the general form of the sums derived from the 
Sine function Feynman rules, we can see that while the 
overall sum is A-dimensional, the summand can typically 
be specified by fewer than N independent parameters. In 
general, the sums take the form 

iAh = E'"EE"'Ep( fci )'"^ ;v ) x 



f(di, ■ ■ ■ ,d p , d p+ x, • • • , d q ). 



(53) 



Rearranging the ki so that {hj+i, • • • , fc/v} only appear 
in a subset, {d Pl , • • • , d q } of {d%, • • • , d q }, we can write 
the innermost N — j sums as a function of {d%, • ■ ■ , d p }, 
F(d±, • • • , dp), where 

F(di,---,d p ) = ^2 •"^2p(kj+i)---p(k N ) X 

f(d 1 ,---,d p ,d p+1 ,---,d q )(54) 



Consequently, our sum is now 

^ = £-£p(*i)-p(Wi 



, d p 



(55) 



This may not seem like an improvement, since we still 
have to evaluate the F(dx, • ■ • , dp). However, by evaluat- 
ing F for a range of values of {dx, • • • , d p } we can construct 
an interpolation table which will yield F to some speci- 
fied accuracy. Constructing the interpolation table will 
typically require exp (c'(7V — j + q — p) operations, while 
the evaluation of the remaining outermost sums will take 
a further exp (c"j) operations. Provided N — j + p — q is 
less than N, the combined process will most likely be more 
efficient than evaluating the N dimensional sum directly. 

We now illustrate this process by the applying it to 
the increasingly familiar sunset diagram. If we focus on 
the renormalized quantity E, we can write equation ( |44| ) 
as 

S„(p) = ^E p(kx)p(k 2 )F(d) (56) 



Figure 10: The three loop diagram that contributes to the 
A0 propagator, with momentum, p. 
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Figure 11: We plot E$j(p) for h = 0.8, m =_j_ and A 2 = 



16. The prefactor of m 2 /(47r) 6 in equation (£ 
removed from these numerical results. 



has been 



where 



F(d) 



P(k 3 ) 



,3 



iMfc 3 )r 

p Z 



exp 



m 2 d+l/c(k 3 ) 



1 



m 2 d + 1/0(^3) 



1 



+ 



c(kx) c(fc 2 ) 



(57) 
(58) 



In practice, by computing F(d) for several hundred values 
of d we calculate to within several parts in 10 10 (with 
h = 0.4) in less than l/5th of the time needed for the 
direct evaluation of the three dimensional sum. 

Another way to speed the evaluation of the diagrams 
would be to exploit the symmetry properties of the sums 
themselves, which reflects the structure of the underlying 
diagram. For instance, for a given {ki,k 2 ,k 3 } the value 
of the general term in is not changed by permuting 
the ki. By using this knowledge it would be possible to 
reduce the number of distinct sets of ki,k 2 , k 3 for which 
the general term in the sum had to be computed. 
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5.3 Higher Order Diagrams 

Until now, we have used the sunset diagram to illustrate 
the Sine function approach to Feynman integrals. How- 
ever, our aim is to develop a method which applies to arbi- 
trary topologies. We now evaluate several specific higher 
order diagrams, although for convenience, we focus on 
contributions to the two-point function. 

These diagrams typically contain sub-divergences, so 
we cannot regularize them using the prescription we em- 
ployed with the sunset diagram. In order to focus on the 
Sine function method itself, in this initial paper we use an 
explicit cut-off to render the higher order diagrams finite, 
and compare the Sine function results to Monte Carlo in- 
tegrations performed with the VEGAS algorithm 00]. 

5.3.1 Case Study II: Third Order 

Consider the following three loop graph, 

£ (3) (zi -x 2 ) = J d 4 yG( Xl - y) 2 G(x 2 - y) 2 G( Xl - x 2 ), 

(59) 

whose topology is illustrated in fig. (JIQ) . Applying the 
Sine function Feynman rules gives: 



m 2 h 5 



p(ki) ■ ■ -p(k 5 ) 



(4tt) 6 ^ (d 1 d 2 + (d 1 +d 2 )/c(k 1 )Y 



exp 



di + d 2 



m 2 d\d 2 + (di + d 2 )/c{ki) 



where 



di 



c(k 2 ) c(fc 3 )' 



d 2 = 



1 



1 



c(fc 4 ) c(fc 5 ) 



, (60) 



(61) 



As with the sunset diagram, h can be taken close to unity 

(3) (3) 

without inducing a large difference between and SW, 

and SX J (p) can be evaluated to six or seven significant 
figures in approximately twenty seconds. Since we have 
no independent estimate of T, A (p) , we heuristically de- 
termine the number of significant figures in our result by 
varying both h and the size of the largest discarded terms 
in the multiple sum. Hopefully, a more careful analysis 
will give an a priori understanding of the truncation error 
in the evaluation of the series, and the h dependence of 
the error in I Ah- 

The results obtained from equation ( pO] ) can be inde- 
pendently checked by using VEGAS to evaluate the five 

(3) 

dimensional integral representation of T, A . As well as 
confirming the Sine function calculation, this also allows 




Figure 12: A four loop diagram that contributes to the 
X(j) 4 propagator, with momentum, p. 

us to compare the time required to evaluate the Sine func- 
tion representation with the Monte Carlo algorithm. We 
find that to obtain four or five significant figures of E(») 
using VEGAS takes around 50 times longer than it took 
us to find seven significant figures of £( 3 ) from I Ah . 

We do not place wish to place too much weight on this 
comparison. In practice, it is very unlikely that any I a 
would be evaluated directly, as it is almost always possible 
to perform some analytical simplifications. On the other 
side of the ledger, we have also not attempted make ana- 
lytical improvements to I Ah- Moreover, VEGAS is a ma- 
ture and well-tested algorithm. Conversely, we are using 
a very simple approach to evaluate the I Ah , arid the pre- 
vious results suggest that a more sophisticated algorithm 
would significantly reduce the time needed to evaluate the 
sums. Finally, the 7a are, in general, cut-off dependent, 
whereas one usually wishes to calculate cut-off indepen- 
dent quantities. These can be combinations of different 
I A and their derivatives, in the limit that A — > oo. This 
process of extracting cut-off independent quantities is pre- 
sumably similar for the Sine function and integral forms, 
but will generally lead to both integrands and summands 
which are more complicated that those which appear in 
the I a and I Ah- 

However, while they are based on an artificial compar- 
ison, these results do suggest that Sine function methods 
may lead to an useful approach for evaluating Feynman 
diagrams, which is intrinsically faster than one which re- 
lies on the use of Monte Carlo integration. Moreover, any 
advantage will become more pronounced when significant 
levels of numerical precision are required. 

5.3.2 Case Study III: Fourth Order 

The final diagram we consider here is the four loop con- 
tribution to the propagator, shown in fig. (|l^), which can 



14 



1.2 ■ 
1 

0.8 
0.6 
0.4 
0.2 



Figure 13: We plot S^O) for h = °- 9 , m=l and A 2 = 
16. As is the case with the 3-loop example we consider, 
this quantity is explicitly cut-off dependent, and we have 
suppressed the m 2 /(47r) 8 prefactor that appears in the 
Sine function representation. 



be represented as the following integral over propagators, 
S (4) (xi-x 2 ) = I d%d i y 2 G(x 1 - x 2 )G(yi ~ y 2 ) 2 x 



G(xi - yi)G(x! - y 2 ) x 

G(x 2 - yi)G{x 2 - y 2 ). (62) 

It is again straightforward to apply the Sine function 
Feynman rules and take the Fourier transform to find 
^aa(p)' an< ^ the results of evaluating it for a specific set 
of A, m and p are shown in fig. (|l3|). 

In this case, VEGAS takes around 20 minutes to find 
a result to within 1 part in 10 4 , whereas the Sine func- 
tion method returns a value accurate to 1 part in 10 6 
in approximately five minutes. Again, the Sine function 
approach outperforms VEGAS by around two orders of 
magnitude, although the seven-dimensional sum is ap- 
preciably more expensive to evaluate than the previous 
five-dimensional example. 

6 Conclusions and Discussion 

In this paper we have developed the Sine function repre- 
sentation, which associates arbitrary Feynman integrals 
with multi-dimensional, infinite sums. The basis of this 
representation is an approximation to the propagator, de- 
rived using the theory of generalized Sine functions. We 
have used the Sine function representation to develop a 
new approach to numerically evaluating Feynman inte- 
grals. This method is simple to implement and is po- 



tentially faster than Monte Carlo algorithms, which are 
currently used to evaluate most Feynman integrals that 
cannot be performed analytically. 

We presented three "case studies" which illustrate the 
properties of the Sine function representation. Firstly, 
we evaluated the two loop sunset diagram, and compared 
the Sine function representation to explicit evaluations of 
the corresponding integral. We showed that we can easily 
compute results that are accurate to within a few parts 
in 10 10 , and that for diagrams with no sub-divergences it 
is straightforward to extract regularization independent 
quantities from the Sine function representation. We then 
evaluated representative third and fourth order diagrams 
from the propagator expansion of X(j) 4 field theory. 

For the three and four loop diagrams we evaluated, 
we contrasted the Sine function representation with direct 
calculations of the corresponding integrals with VEGAS, 
showing that the Sine function representation is signif- 
icantly more efficient than VEGAS. This advantage in- 
creases as the desired level of accuracy is raised, and can 
easily amount to several orders of magnitude in execution 
time. However, this direct comparison between Monte 
Carlo and Sine function methods is somewhat artificial 
since the quantities computed are cut-off dependent, and 
we did not attempt to analytically simplify the integrals 
(or the sums) before evaluating them. Consequently, it re- 
mains to be seen whether the Sine function representation 
will speed up realistic calculations. 

In order to test the Sine function representation in 
practical situations, we plan to investigate a number of 
topics. We have already used the Sine function represen- 
tation to evaluate the three-loop master diagrams whose 
analytical properties are discussed in detail by Broad- 
hurst pl| . In the course of this work |L6| we explicitly ex- 
tended the Sine function representation to diagrams with 
massless lines, and reproduced the analytical values of 
the master diagrams to better than 1 part in 10 10 . In ad- 
dition to describing a renormalization procedure for the 
Sine function representation of arbitrary diagrams, our 
immediate priority is to include fermion and vector fields, 
and thus to evaluate diagrams from QED and electroweak 
theory. We have calculated several simple QED diagrams 
with Sine function techniques fl7|| , and we are currently 
working on the systematic application of the Sine function 
representation to arbitrary QED diagrams. 

We have used a straightforward procedure to evaluate 
the multi-dimensional sums derived from the Sine func- 
tion Feynman rules. In Section 4 we discussed a vari- 
ety of approaches to speeding up the computation of the 
sums, and it is possible that the numerical performance of 
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the Sine function representation can be significantly im- 
proved. Moreover, we have not yet attempted to analyti- 
cally simplify Sine function representations before evalu- 
ating then numerically. Consequently, the CPU time we 
currently need to evaluate the sum only provides an upper 
limit on what can be achieved. 

Obtaining the Sine function representation for a spe- 
cific topology amounts to evaluating one or more Gaus- 
sian integrals, and the process is easily automated with 
any computer algebra package. If the Sine function repre- 
sentation is useful in realistic situations, we expect it will 
be straightforward to integrate it with existing tools for 
automatically evaluating diagrams. Moreover, applying 
the Sine function Feynman rules to diagrams which con- 
tain several fields with different masses is no more com- 
plicated than evaluating the same topology with identical 
lines. Diagrams with more than two external lines also 
present no new difficulties, although in this case more 
than one Fourier transform is needed to move from coor- 
dinate to momentum space, since an n-point diagram has 
n — 1 independent external momenta. 

Strictly speaking the Sine function representation is 
not just a "better way to do integrals" , since when evalu- 
ated to arbitrary precision, the Sine function representa- 
tion converges to an approximation to the integral, rather 
than the integral itself. However, this approximation can 
be made arbitrarily accurate, so in practice evaluating 
the Sine function representation is equivalent to a direct 
computation of the integral. Moreover, several other ap- 
proaches lead to expressions for Feynman integrals that 
are written as infinite sums p8|-pl[|. However, these meth- 
ods yield expressions for the exact integrals, and are there- 
fore not equivalent to the Sine function representation. 

While many approaches to numerical integration out- 
perform Monte Carlo techniques in one or two dimensions, 
their efficiency scales typically very poorly with the di- 
mension of the integral. For the integrals derived from 
Feynman diagrams, the Sine function representation pro- 
vides an exception to this general rule. For example, the 
4th order diagram we studied can be expressed as a seven 
dimensional integral, but a Monte Carlo integration with 
VEGAS takes much longer than the Sine function repre- 
sentation to achieve the same level of accuracy. The task 
before us now is to harness the mathematical properties 
of the Sine function representation, and apply them to 
solving physical problems. 



A Numerical Details 

The numerical results described in this paper we all ob- 
tained on the same Sun workstation, with 250MHz Ul- 
trasparc II CPUs. We refrain from giving precise timings, 
since these depend strongly on the specific combination of 
hardware and software, and the timings we do give reflect 
the use of aggressive compiler optimization settings. The 
codes are implemented in Fortran 77. Sample codes (in- 
cluding those used to perform the calculations described 

in this paper) are available at the following URL: 

http : //www .net .brown . edu/people/easther /feynman 
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